ASMBI-Paper about POMDP
Enhao Liu
June 11, 2019
1 Library Settings
Note: we must prevent the scientific notation issue. For example, when using paste("abc", 10000), we want the results to be "abc 10000" instead of "abc 1e+05". Hence, using options("scipen" =10) would deal with this issue. The reason of doing this is that the solve_POMDP() would first call write_POMDP() which convert the transition probabilities, observation probabilities and rewards to be a file.POMDP. The key of method is using paste() function. If there is no prevention for scientific notation, then the results from paste() would be some characters containing "1e-09", thereby causing that the Initilization POMDP part have several syntax error.
2 Introduction
In this R markdown, we will redo the development of cyber vulnerabilities policies based on POMDP, which is in our submitted ASMBI paper. The parameters of POMDP have been estimated and shown in the paper.
3 POMDP model
3.1 Discount, States, Actions, Observations
discount<-0.99
states<-c("low" , "medium" , "high", "critical")
actions<-c("auto-patching" , "research-accept", "research-compensate", "remediation")
observations<-c("none", "low", "medium","high","critical")
# states <- c("s1", "s2", "s3", "s4")
# actions <- c("a1", "a2", "a3", "a4")
# observations <- c("obs1", "obs2", "obs3", "obs4", "obs5")3.2 Transitions Prob., Observations Prob. and Rewards
3.2.1 Load
First, we load the estimated parameters.
#Transition for OSs
T_linux<-read.csv("Data for POMDP/Parameters_transition_linux.csv")
T_windows<-read.csv("Data for POMDP/Parameters_transition_windows.csv")
T_other<-read.csv("Data for POMDP/Parameters_transition_other.csv")
#Obervation
Prop_disappeared<-read.csv("Data for POMDP/Disapeared Observation Proportion_approach2.csv")
#Reward
R_all<-read.csv("Data for POMDP/Parameters_reward.csv")3.2.2 Data structure for POMDP solver
Next, we modify the data structure for POMDP solver.
(1) Transition Probabilities
transition_prob is a \(OS\times A\times S\times S'\) tensor.
n_state<-length(states) # number of states
n_obs<-length(observations) # number of observations
n_action<-length(actions) # number of actions
#convert transition prob dataframe to be matrix
T_linux<-T_linux%>%as.matrix()%>%unname()
T_windows<-T_windows%>%as.matrix()%>%unname()
T_other<-T_other%>%as.matrix()%>%unname()
#construct transition prob for POMDP: list data structure
transition_prob<-list(T_linux, T_windows, T_other)
names(transition_prob)<-c("Linux", "Windows", "Other")
transition_prob<-lapply(transition_prob, function(x){
prob<-lapply(1:n_action, function(i){
x[ ,(1+(i-1)*n_state):(i*n_state)]
})
names(prob)<-actions
return(prob)
})
# example of transition prob: Linux
transition_prob$Linux## $`auto-patching`
## [,1] [,2] [,3] [,4]
## [1,] 0.9171 0.0800 0.0020 0.0009
## [2,] 0.0031 0.9905 0.0059 0.0005
## [3,] 0.0000 0.0032 0.9905 0.0063
## [4,] 0.0000 0.0000 0.0032 0.9968
##
## $`research-accept`
## [,1] [,2] [,3] [,4]
## [1,] 0.9963 0.0037 0.0000 0.0000
## [2,] 0.4299 0.5664 0.0037 0.0000
## [3,] 0.0316 0.3983 0.5664 0.0037
## [4,] 0.0091 0.0604 0.0205 0.9100
##
## $`research-compensate`
## [,1] [,2] [,3] [,4]
## [1,] 1.0000 0.0000 0.0000 0
## [2,] 1.0000 0.0000 0.0000 0
## [3,] 0.0315 0.9648 0.0037 0
## [4,] 0.0091 0.0604 0.9305 0
##
## $remediation
## [,1] [,2] [,3] [,4]
## [1,] 0.5 0.5 0 0
## [2,] 0.5 0.5 0 0
## [3,] 0.5 0.5 0 0
## [4,] 0.5 0.5 0 0(2) Observation Probabilities
observation_prob is a \(OS\times Criticality\times A\times S'\times O\).
# modify the Critical column to be categorical
Prop_disappeared2<-Prop_disappeared%>%
mutate(Critical = paste("crit",Critical,sep = "_"))
Prop_disappeared2#construct observation prob for POMDP: list data structure
Prop_disappeared_list<-plyr::dlply(Prop_disappeared2, "OS.Type")
observation_prob<-lapply(Prop_disappeared_list, function(x){
x<-plyr::dlply(x, "Critical")
x<-lapply(x, function(y){
prop_missing<-y%>%select(Average)%>%pull()
obs_prob<-matrix(data = 0, nrow = n_state, ncol = n_obs)
obs_prob[,1]<-prop_missing #the first column = proportion of missing
for (j in 2:n_obs){ #the second and after = 1-proportion of missing
obs_prob[j-1,j] = 1-prop_missing
}
#assign to each action
prob<-lapply(1:n_action, function(i){
obs_prob
})
names(prob)<-actions
return(prob)
})
})
#Example of observation prob: Linux and critical=False
observation_prob$Linux$crit_FALSE## $`auto-patching`
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.3582727 0.6417273 0.0000000 0.0000000 0.0000000
## [2,] 0.3582727 0.0000000 0.6417273 0.0000000 0.0000000
## [3,] 0.3582727 0.0000000 0.0000000 0.6417273 0.0000000
## [4,] 0.3582727 0.0000000 0.0000000 0.0000000 0.6417273
##
## $`research-accept`
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.3582727 0.6417273 0.0000000 0.0000000 0.0000000
## [2,] 0.3582727 0.0000000 0.6417273 0.0000000 0.0000000
## [3,] 0.3582727 0.0000000 0.0000000 0.6417273 0.0000000
## [4,] 0.3582727 0.0000000 0.0000000 0.0000000 0.6417273
##
## $`research-compensate`
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.3582727 0.6417273 0.0000000 0.0000000 0.0000000
## [2,] 0.3582727 0.0000000 0.6417273 0.0000000 0.0000000
## [3,] 0.3582727 0.0000000 0.0000000 0.6417273 0.0000000
## [4,] 0.3582727 0.0000000 0.0000000 0.0000000 0.6417273
##
## $remediation
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.3582727 0.6417273 0.0000000 0.0000000 0.0000000
## [2,] 0.3582727 0.0000000 0.6417273 0.0000000 0.0000000
## [3,] 0.3582727 0.0000000 0.0000000 0.6417273 0.0000000
## [4,] 0.3582727 0.0000000 0.0000000 0.0000000 0.6417273Here, we consider the error scenarios to vary observation matrix. For example, let error be 10%.
- For actual \(s' = Low\), observing medium, high or critical state is impossible. So such error can not affect this.
- For actual \(s' = Medium\), we might observe \(Low\) with 10% probability. Then,
\[ P_{error}(o = Medium|s' = Medium) = P(o = Medium|s'= Medium)*(1-10\%) \\ P_{error}(o = Low|s' = Medium) = P(o = Medium|s'= Medium)*(10\%) \]
- For actual \(s' = High\), we might observe \(Low\) or \(Medium\) with 10% probability. Then,
\[ P_{error}(o = High|s' = High) = P(o = High|s'= High)*(1-10\%) \\ P_{error}(o = Medium|s' = High) = P(o = High|s'= High)*(10\%/2)\\ P_{error}(o = Low|s' = High) = P(o = High|s'= High)*(10\%/2) \]
- For actual \(s' = High\), we might observe \(Low\) or \(Medium\) with 10% probability. Then,
\[ P_{error}(o = High|s' = High) = P(o = High|s'= High)*(1-10\%) \\ P_{error}(o = Medium|s' = High) = P(o = High|s'= High)*(10\%/2)\\ P_{error}(o = Low|s' = High) = P(o = High|s'= High)*(10\%/2) \]
- For actual \(s' = Critical\), we might observe \(Low\), \(Medium\) or \(Critical\) with 10% probability. Then,
\[ P_{error}(o = Critical|s' = Critical) = P(o = Critical|s'= Critical)*(1-10\%) \\ P_{error}(o = High|s' = Critical) = P(o = Critical|s'= Critical)*(10\%/3) \\ P_{error}(o = Medium|s' = Critical) = P(o = Critical|s'= Critical)*(10\%/3)\\ P_{error}(o = Low|s' = Critical) = P(o = Critical|s'= Critical)*(10\%/3) \]
The above procedure has been put into the following function Consider_Error_Observation
# Consider error in observation matrix
Consider_Error_Observation <- function(obs_matrix, error_scan = 0){
obs_matrix <- lapply(obs_matrix, function(mat){
for (i in 2:nrow(mat)){
if(i == 2){
dummy_value <- mat[i, i+1]
mat[i, i] <- dummy_value*error_scan
mat[i, i+1] <- dummy_value*(1-error_scan)
}
if (i==3){
dummy_value <- mat[i, i+1]
mat[i,i-1] <- dummy_value*(error_scan/2)
mat[i,i] <- dummy_value*(error_scan/2)
mat[i, i+1] <- dummy_value*(1-error_scan)
}
if (i==4){
dummy_value <- mat[i, i+1]
mat[i,i-2] <- dummy_value*(error_scan/3)
mat[i,i-1] <- dummy_value*(error_scan/3)
mat[i,i] <- dummy_value*(error_scan/3)
mat[i, i+1] <- dummy_value*(1-error_scan)
}
}
return(mat)
})
return(obs_matrix)
}
# Consider_Error_Observation(obs_matrix = observation_prob$Linux$crit_FALSE,
# error_scan = 0.1)(3) Reward
reward is a \(OS\times Criticality\times Admin \times Dataframe\) where \(Dataframe\) contains columns: action, start-state, end-state, observation and reward.
reward<-R_all%>%mutate(start.state = rep(states, times = nrow(R_all)/n_state),
end.state = "*",
observation = "*")%>%
gather(key = "action", value = "reward", -OS,-Critical,-Admin,-start.state,-end.state,-observation)%>%
select(OS, Critical, Admin, action, start.state, end.state,observation, reward)
#match the action names
actions_old<-reward%>%distinct(action)%>%pull()
reward<-reward%>%mutate(action = plyr::mapvalues(action, from = actions_old, to = actions))
#construct a list structure:
#Step1 : convert the logical value in Critical and Admin columns to be character values
reward<-reward%>%
mutate(Critical = paste("crit",Critical,sep = "_"))%>%
mutate(Admin = paste("admin", Admin, sep = "_"))
#Step2: split by OS to make a list
reward_list<-plyr::dlply(reward, "OS")
#Step3: split the OS list by Critical and then by Admin. After that, removing the OS, Critical and Admin columns.
reward_list<-lapply(reward_list, function(x){
x<-plyr::dlply(x, "Critical") #split each OS list by "Critical or not"
x<-lapply(x,function(y){
plyr::dlply(y, "Admin") #split nested OS and Critical list by Admin
})
x<-lapply(x, function(y){ #remove OS,Critical and Admin columns in nested OS,Critical and Admin dataframe
lapply(y, function(z){
z%>%select(-OS, -Critical, -Admin)%>%
mutate_if(is.character, as.factor)
})
})
return(x)
})
#Example of reward: Linux, critical=False, admin = False
reward_list$Linux$crit_FALSE$admin_FALSE4 Solving POMDPs Functions
There are some issues in solve_POMDP() from pomdp package. Specifically, the issues happen in $pg. Therefore, we have to modify the solve_POMDP() function. In addition, some functions have been made.
4.1 find_pomdpsolve function
4.2 belief_update function
Every time agent taking an action \(a\) and observes \(o\), its belief is updated by Bayes rule:
\[ b^{ao}(s') = \frac{p(o|s',a)}{p(o|b,a)}\sum\limits_{s\in S}p(s'|s,a)*b(s) \]
where \(p(s'|s,a)\) and \(p(o|s',a)\) are defined by the model parameters transition probability and observation probability, and
\[ p(o|b,a) = \sum\limits_{s'\in S}p(o|s', a)\sum\limits_{s\in S}p(s'|s,a)*b(s) \]
whichh is normalization term.
Here, we define belief_update() function to implement the above formulation.
belief_update <- function(b_0, model, act_index, obs_index){
if (class(model)!= "POMDP_model"){
stop("the 'model' argument must be constructed by POMDP() funtion")
}else {
n_observations <- length(model$observations) # number of observations
n_state <- length(model$states) # number of states
# if (act_index > length(model$actions)){
# stop("Error in act_index: it is out of bound")
# }
# if(obs_index > n_observations){
# stop("Error in obs_index: it is out of bound")
# }
transition_prob_act <- model$transition_prob[[act_index]]%>%as.matrix() # transition prob of action
observation_prob_act <- model$observation_prob[[act_index]]%>%as.matrix() #observation prob of action
if (observation_prob_act[1,1] == "uniform"){
observation_prob_act <- matrix(1/n_observations, nrow = n_state, ncol = n_observations)
}
# compute the denominator (normalization term)
sum_1<-0
sum_2<-0
for (s_end in 1:n_state){
for (s_start in 1:n_state){
x<-b_0[s_start]*transition_prob_act[s_start, s_end]
sum_1<-sum_1+x
}
y<-sum_1*observation_prob_act[s_end, obs_index]
sum_2<-sum_2+y
#print(sum_1)#debug
}
p_normalizer<-sum_2
# update belief
if(p_normalizer == 0){
b_update <- rep(NA_real_, n_state)
cat(paste("Warnings in ", model$name,": The belief update is NAs because it is impossible to observe",
model$observations[obs_index],
"observation under action",
model$actions[act_index],
"when the start belief is (",
paste(b_0, collapse = ","), ")",
sep = " "))
return(b_update)
}else{
b_update<-b_0
sum<-0
for (s_end in 1:n_state){
for (s_start in 1:n_state){
x<-b_0[s_start]*transition_prob_act[s_start, s_end]
sum<-sum+x
}
b_update[s_end]<-observation_prob_act[s_end, obs_index]*sum/p_normalizer
}
return(b_update)
}
}
}
#belief_update(b_0 = c(1,0,0,0), model = model_POMDPs[[1+12]], act_index = 1, obs_index = 3)4.3 create_PG function
Since the pg dataframe has some problem, here we are using create_PG function to re-construct the policy graph.
The input of this function is:
model: it is POMDP model which has been constructed byPOMDP()functionalpha: it is a dataframe that contains alpha vectors. Each row refers to an alpha vector. The number of columns is equal to the number of states.action_index_alpha: it is a vector of actions index. Each action is associated with the row of alpha dataframe.belief_matrix: it is a matrix. Each row refers to a vector of belief, for example, (0.25, 0.25, 0.25, 0.25). The number of columns is equal to the number of states.
create_PG <- function(model, alpha, action_index_alpha, belief_matrix){
if (class(model)!= "POMDP_model"){
stop("Error in model: the 'model' argument must be constructed by POMDP() funtion")
}else {
if(is.null(belief_matrix)){
pg <- NULL
return(pg)
}else{
alpha_matrix <- as.matrix(alpha)
n_alpha <-dim(alpha_matrix)[1]
if(length(action_index_alpha)!=n_alpha){
stop("Error in action_index_alpha: the length of this vector has to be the same as the number of rows of alpha")
}
# PG construction:
# initialize optimal action index for each row of belief matrix
act_opt_index <- rep(0, dim(belief_matrix)[1])
# initialize PG
n_observations <- length(model$observations)
PG_int <- matrix(0, nrow = dim(belief_matrix)[1], ncol = 2 + n_observations)
# for each row of belief matrix, find the index of alpha vector that gives maximum value
# and then find the opt action index
for (i in 1:dim(belief_matrix)[1]){
if (is.na(belief_matrix[i, 1])){ #if the belief points are NAs
PG_int[i, 1] <-NA_integer_
}else{
alpha_opt_index <- which.max(alpha_matrix %*% belief_matrix[i,])
act_opt_index[i] <- action_index_alpha[alpha_opt_index]
PG_int[i,1] <- alpha_opt_index
}
}
PG_int[ ,2] <- act_opt_index
# for each row of belief matrix, get the alpha opt index for each observation
for (i in 1:dim(belief_matrix)[1]){
for (j in 1:n_observations){
if (is.na(belief_matrix[i, 1])){ #if the belief points are NAs
PG_int[i, j+2] <- NA_real_
}else{
b_update <- belief_update(b_0 = belief_matrix[i, ],
model = model,
act_index = act_opt_index[i],
obs_index = j)
if (is.na(b_update[1])){ #if the updated belief is NA
PG_int[i, j+2] <- NA_real_
}else{
PG_int[i, j+2] <- which.max(alpha_matrix %*% b_update)
}
}
}
}
#rename and construct PG as dataframe
pg <- as.data.frame(PG_int)
# renamig the columns
colnames(pg)[1] <- "belief_state"
colnames(pg)[2] <- "action"
for (i in 1:n_observations) {
colnames(pg)[i+2] <- model$observations[i]
}
# renaming the actions
for (i in 1:length(model$actions)) {
pg[pg[,2]==i,2] <- model$actions[i]
}
return(pg)
}
}
}
# create_PG(model = LinuxProblem2,
# alpha = Linux_sol$alpha,
# action_index_alpha = action_index_alpha,
# belief_matrix = belief_prop_matrix)4.4 create_PG2 function
This one will generate PG with the action name or action index
create_PG2 <- function(model, alpha, action_index_alpha, belief_matrix,
get_action_name = TRUE){
if (class(model)!= "POMDP_model"){
stop("Error in model: the 'model' argument must be constructed by POMDP() funtion")
}else {
if(is.null(belief_matrix)){
pg <- NULL
return(pg)
}else{
alpha_matrix <- as.matrix(alpha)
n_alpha <-dim(alpha_matrix)[1]
if(length(action_index_alpha)!=n_alpha){
stop("Error in action_index_alpha: the length of this vector has to be the same as the number of rows of alpha")
}
# PG construction:
# initialize optimal action index for each row of belief matrix
act_opt_index <- rep(0, dim(belief_matrix)[1])
# initialize PG
n_observations <- length(model$observations)
PG_int <- matrix(0, nrow = dim(belief_matrix)[1], ncol = 2 + n_observations)
# for each row of belief matrix, find the index of alpha vector that gives maximum value
# and then find the opt action index
for (i in 1:dim(belief_matrix)[1]){
if (is.na(belief_matrix[i, 1])){ #if the belief points are NAs
PG_int[i, 1] <-NA_integer_
}else{
alpha_opt_index <- which.max(alpha_matrix %*% belief_matrix[i,])
act_opt_index[i] <- action_index_alpha[alpha_opt_index]
PG_int[i,1] <- i-1
}
}
if (get_action_name){
PG_int[ ,2] <- model$actions[act_opt_index]
}else{
PG_int[ ,2] <- act_opt_index
}
# for each row of belief matrix, get the alpha opt index for each observation
for (i in 1:dim(belief_matrix)[1]){
for (j in 1:n_observations){
if (is.na(belief_matrix[i, 1])){ #if the belief points are NAs
if(get_action_name){
PG_int[i, j+2] <- NA_character_
}else{
PG_int[i, j+2] <- NA_integer_
}
}else{
b_update <- belief_update(b_0 = belief_matrix[i, ],
model = model,
act_index = act_opt_index[i],
obs_index = j)
if (is.na(b_update[1])){ #if the updated belief is NA
if(get_action_name){
PG_int[i, j+2] <- NA_character_
}else{
PG_int[i, j+2] <- NA_integer_
}
}else{
alpha_opt_index <- which.max(alpha_matrix %*% b_update)
act_index <- action_index_alpha[alpha_opt_index] # get action index
if (get_action_name){
PG_int[i, j+2] <- model$actions[act_index]
}else{
PG_int[i, j+2] <- act_index
}
}
}
}
}
#rename and construct PG as dataframe
pg <- as.data.frame(PG_int)
# renamig the columns
colnames(pg)[1] <- "belief_scenario"
colnames(pg)[2] <- "action"
for (i in 1:n_observations) {
colnames(pg)[i+2] <- model$observations[i]
}
# # renaming the actions
# for (i in 1:length(model$actions)) {
# pg[pg[,2]==paste(i),2] <- model$actions[i]
# }
return(pg)
}
}
}
# create_PG2(model = model_POMDPs_solved[[5]]$model,
# alpha = model_POMDPs_solved[[5]]$solution$alpha,
# action_index_alpha = model_POMDPs_solved[[5]]$solution$action_index_alpha,
# belief_matrix = matrix(c(0.25, 0.25, 0.25, 0.25,
# 1, 0, 0, 0,
# 0, 1, 0, 0,
# 0, 0, 1, 0,
# 0, 0, 0, 1), nrow = 5, ncol = 4, byrow = TRUE),
# get_action_name = TRUE)4.5 solve_POMDP2 function
This function is modified based on solve_POMDP. The key difference is that we re-calculate the pg dataframe based on belief_update function, create_PG function, and alpha vector instead of reading the .pg file from the results of pomdp-solve.exe.
solve_POMDP2 <- function(
model,
horizon = NULL,
method = "grid",
parameter = NULL,
verbose = FALSE) {
### DEBUG
#verbose <- TRUE
###
methods <- c("grid", "enum", "twopass", "witness", "incprune")
### Not implemented: "linsup", "mcgs"
method <- match.arg(method, methods)
### write model to file
file_prefix <- tempfile(pattern = "model")
pomdp_filename <- paste0(file_prefix, ".POMDP")
pomdp::write_POMDP(model, pomdp_filename)
### running the POMDP code
if(!is.null(parameter)) {
paras <- sapply(names(parameter), FUN = function(n) paste0("-", n, " ", parameter[[n]]))
} else paras <- ""
solver_output <- system2(find_pomdpsolve(), #invoke a system command
args = c(paste("-pomdp", pomdp_filename),
paste("-method", method),
(if(!is.null(horizon)) paste("-horizon", horizon) else ""),
paras,
"-fg_save true"),
stdout = TRUE,
stderr = TRUE,
wait = TRUE
)
## the verbose mode: printing all the outputs from pomdp solver
if(verbose) cat(paste(solver_output, "\n"))
### importing the outputs and results (pomdp-solve adds the PID to the file prefix)
file_prefix <- gsub("^o = (.*)\\s*$", "\\1", solver_output[16])
## Creating result files' names and extensions
pomdp_filename <- paste0(file_prefix, ".POMDP")
pg_filename <- paste0(file_prefix, ".pg")
belief_filename <- paste0(file_prefix, ".belief")
alpha_filename <- paste0(file_prefix, ".alpha")
## importing alpha file
number_of_states <- length(model$states)
number_of_observations <- length(model$observations)
number_of_actions <- length(model$actions)
observations <- model$observations
actions <- model$actions
states <- model$states
start <- model$start
alpha <- read.table(alpha_filename, header = FALSE, sep = "\n")
alpha <- as.matrix(alpha)
toDel <- seq(1,dim(alpha)[1],2) #the row index which contains the value of action
action_index_alpha <- alpha[toDel, 1]%>%as.numeric() # the actions index associated with each alpha vector
action_index_alpha <- action_index_alpha +1
alpha <- alpha[-toDel,]
alpha <- unlist(strsplit(alpha, " "))
alpha<- matrix(as.numeric(alpha), ncol = number_of_states, byrow = TRUE)
alpha_matrix <- alpha
alpha <- as.data.frame(alpha)
# renaming the columns
for (i in 1:number_of_states) {
colnames(alpha)[i] <- paste("coeffecient", i)
}
## importing pg file
# pg <- read.table(pg_filename, header = FALSE, sep = " ", quote = "\"",
# dec = ".", na.strings = NA, numerals = "no.loss")
# pg <- as.matrix(pg)
# pg <- pg[,c(-3,-dim(pg)[2])] #there 2 NA columns that need to be removed
# pg <- pg+1 #index has to start from 1 not 0
# pg_matrix <- pg
# pg <- as.data.frame(pg)
# if (dim(pg)[2]==1 ) {
# pg <- t(pg)
# pg <-as.data.frame(pg)
# }
## importing belief file if it exists
if(file.exists(belief_filename)) {
belief <- read.table(belief_filename)
belief_matrix <- as.matrix(belief)
belief <- as.data.frame(belief_matrix)
# renaming the columns
for (i in 1:number_of_states) {
colnames(belief)[i] <- states[i]
}
## finding the respective proportions for each line (node)
for (i in 1:dim(belief_matrix)[1]) {
belief[i,number_of_states+1]<- which.max(alpha_matrix %*% belief_matrix[i,])
}
colnames(belief)[number_of_states+1] <- "belief_state"
belief_proportions <- alpha-alpha
colnames(belief_proportions) <- colnames(belief)[1:number_of_states]
for (i in 1:dim(alpha_matrix)[1]) {
c <- 0
for (j in 1: dim(belief_matrix)[1]) {
if (belief[j,ncol(belief)]==i) {
c <- c + 1
belief_proportions[i,] <- belief_proportions[i,]+belief_matrix[j,]
}
}
belief_proportions[i,] <- belief_proportions[i,]/c
}
} else {
belief <- NULL
belief_proportions <- NULL
}
## ============Construct PG=====================
if (is.null(belief_proportions)){
pg <- NULL
}else{
belief_prop_matrix <- as.matrix(belief_proportions)
pg <- create_PG(model = model,
alpha = alpha,
action_index_alpha = action_index_alpha,
belief_matrix = belief_prop_matrix)
}
### outputs and results
## producing the starting belief vector
if (!is.character(start)) {
if (sum(start)==1) {
start_belief <- start
}
}
# if the starting beliefs are given by a uniform distribution over all states
if (sum(start== "uniform")==1) {
start_belief <- rep(1/number_of_states,number_of_states)
} else if (start[1]!="-") { # if the starting beliefs include a specific subset of states
# if the starting beliefs are given by a uniform distribution over a subset of states (using their names)
if (!is.na(sum(match(start , states)))) {
start_belief <- rep(0,number_of_states)
start_belief[match(start,states)] <- 1/length(start)
}
# if the starting beliefs are given by a uniform distribution over a subset of states (using their numbers)
if (!is.character(start)) {
if (sum(start)>=1 & length(start)<number_of_states) {
start_belief <- rep(0,number_of_states)
start_belief[start] <- 1/length(start)
}
}
} else if (start[1]=="-") { # if the starting beliefs exclude a specific subset of states
start_belief <- rep(1/(number_of_states-length(start)+1),number_of_states)
if (is.na(as.numeric(start[2]))) {
start_belief[match(start,states)] <- 0
}
if (!is.na(as.numeric(start[2]))) {
start_belief[start] <- 0
}
}
## calculating the total expected reward for start_belief
initial_belief_state <- which.max(alpha_matrix %*% start_belief)
total_expected_reward <- max(alpha_matrix %*% start_belief)
solution <- structure(list(
method = method,
parameter = parameter,
alpha = alpha,
action_index_alpha = action_index_alpha,
pg = pg,
belief = belief,
belief_proportions = belief_proportions,
total_expected_reward = total_expected_reward,
initial_belief_state = initial_belief_state
), class = "POMDP_solution")
structure(list(model = model,
solution = solution,
solver_output = solver_output
), class = "POMDP")
}
print.POMDP <- function(x, ...) {
cat(
"Solved POMDP model:", x$model$name, "\n",
"\tmethod:", x$solution$method, "\n",
"\tbelief states:", nrow(x$solution$pg), "\n",
paste0("\ttotal expected ", x$model$values, ":"),
x$solution$total_expected_reward,"\n\n"
)
}
print.POMDP_solution <- function(x, ...) {
cat("POMDP solution\n\n")
print(unclass(x))
}5 Simulation of POMDP functions
5.1 Action_map function
# Action_map is used to get the action index for belief under policy
Action_map <- function(belief, policy = "optimal", alpha, action_index_alpha, obs){
# if (policy == "fix"){
# if (belief[3] > 0.5 | belief[4] > 0.5){
# action_index <- 2
# }else if(belief[1] > 0.5 | belief[2] > 0.5){
# action_index <- 1
# }else{
# action_index <- 5 # action5==ignore
# }
# }
if (policy == "fix"){
state_index <- which.max(belief)
if (state_index == 3 | state_index == 4){
action_index <- 2
}else{
action_index <- 1
}
}
# if (policy == "fix"){
# if (belief[3] == 1 | belief[4] == 1){
# action_index <- 2
# }else{
# action_index <- 1
# }
# }
# =====fix policy based on observations====
if (policy == "OSOM"){
if (obs == 4 | obs == 5){
action_index <- 2
}else if(obs == 2 | obs == 3){
action_index <- 1
}else{ #if obs=1, observe nothing, do action5==ignore
action_index <- 5
}
}
# if (policy == "fix_1111"){
# if (obs == 1){
# action_index <- 5
# }else{
# action_index <- 1
# }
# }
if (policy == "optimal"){
alpha_matrix <- as.matrix(alpha)
get_index <- which.max(alpha_matrix %*% belief)
action_index <- action_index_alpha[get_index]
}
if (policy == "random"){
action_index <- sample(1:4, size = 1, replace = TRUE)
}
if (policy == "random_12"){
action_index <- sample(1:2, size = 1, replace = TRUE)
}
if (policy == "e_greedy"){
alpha_matrix <- as.matrix(alpha)
get_index <- which.max(alpha_matrix %*% belief)
action_index_opt <- action_index_alpha[get_index] #optimal action index
action <- c(1:4)
action_exclude <- action[action!= action_index_opt] #remove the optimal action index
action_index_greedy <- sample(action_exclude, size = 1, replace = TRUE)
action_index <- sample(c(action_index_opt, action_index_greedy),
size = 1,
replace = TRUE,
prob = c(0.95, 0.05))
}
return(action_index)
}
# Action_map(belief = c(0,0,0,1),
# policy = "fix",
# alpha,
# action_index_alpha)5.2 Simulate_policy function
The main idea in this Simulate_policy function is to randomly generate \(state, observation\) based on their probabilities under \(action, belief\) at each time step.
Simulate_policy <- function(model_solved, n_trial = 1000, n_horizon = 60,
policy = "fix", b_0 = "random", seed = 123){
states <- model_solved$model$states
actions <- model_solved$model$actions
observations <- model_solved$model$observations
gamma <- model_solved$model$discount
n_state <- length(states)
n_obs <- length(observations)
# policy
if (policy == "optimal"){
alpha <- model_solved$solution$alpha
action_index_alpha <- model_solved$solution$action_index_alpha
}
if (policy == "e_greedy"){
alpha <- model_solved$solution$alpha
action_index_alpha <- model_solved$solution$action_index_alpha
}
# reward matrix
reward <- model_solved$model$reward
reward <- reward %>%
mutate(start.state = fct_relevel(start.state, states))%>%
spread(key = "action", value = "reward")%>%
select(actions)
# observations matrix
observation_prob <- model_solved$model$observation_prob
# transition matrix
transition_prob <- model_solved$model$transition_prob
# # belief matrix
# belief_matrix <- model_solved$solution$belief[, 1:n_state]%>%as.matrix()%>%unname()
# if policy is fixed, add action5 = ignore
if(policy == "fix" | policy == "OSOM"){
#transition_prob[["ignore"]] <- diag(x= 1, nrow = n_state, ncol = n_state)
# transition_prob[["ignore"]] <- matrix(data = c(1, 0, 0, 0,
# 0, 1/3, 1/3, 1/3,
# 0, 0, 0.5, 0.5,
# 0, 0, 0, 1),
# nrow = n_state, ncol = n_state, byrow = TRUE)
transition_prob[["ignore"]] <- matrix(data = 1/n_state,
nrow = n_state, ncol = n_state, byrow = TRUE)
# observation_prob[["ignore"]] <- matrix(data = c(1,0,0,0,0,
# 1,0,0,0,0,
# 1,0,0,0,0,
# 1,0,0,0,0),
# nrow = n_state, ncol = n_obs, byrow = T)
observation_prob[["ignore"]] <- model_solved$model$observation_prob[[1]]
reward <- reward%>%dplyr::mutate(ignore = .[,1]) #action5=ignore cost == action 1 cost
model_solved$model$observation_prob <- observation_prob
model_solved$model$transition_prob <- transition_prob
}
#
# ==========Simulation===============
# run n_trial to get expected rewards
# store results
r_total_trials <- rep(0, n_trial)
action_trials <- matrix(0, nrow = n_trial, ncol = n_horizon)
state_trials <- matrix(0, nrow = n_trial, ncol = n_horizon)
obs_trials <- matrix(0, nrow = n_trial, ncol = n_horizon)
r_trials <- matrix(0, nrow = n_trial, ncol = n_horizon)
set.seed(seed) # for reproduce
#======outer loop: trials start=======
for (t in 1:n_trial){
discount <- 1
r_total <- 0
#====initial version 1:
if (b_0 == "uniform"){
obs <- 1
b <- rep(1/n_state, n_state)
s_start <- sample(c(1:n_state), size = 1, replace = TRUE, prob = b) #random generate state: s_start
}
#===initial version 2: generate b_0 based on random observation at begining
if (b_0 == "random"){
obs <- sample(c(1:n_obs), size = 1, replace = TRUE)
if (obs == 1){ #if observe nothing, randomly select b_0 from belief matrix
# b_index <- sample(c(1:nrow(belief_matrix)), size = 1, replace = TRUE)
# b <- belief_matrix[b_index, ]
b <- rep(1/n_state, n_state)
}else if(obs == 2){
b <- c(1, 0, 0, 0)
}else if(obs == 3){
b <- c(0, 1, 0, 0)
}else if (obs == 4){
b <- c(0, 0, 1, 0)
}else{
b <- c(0, 0, 0, 1)
}
s_start <- sample(c(1:n_state), size = 1, replace = TRUE, prob = b) #random generate state: s_start
}
#====inner loop: periods start======
for (h in 1:n_horizon){
# get action index: a
a <- Action_map(belief = b, policy = policy, alpha, action_index_alpha, obs)
action_trials[t, h] <- a
state_trials[t, h] <- s_start
# random generate next state: s_end
s_end <- sample(c(1:n_state), size = 1, replace = TRUE, prob = transition_prob[[a]][s_start, ])
# random generate observation: obs
obs_trials[t, h] <- obs
obs <- sample(c(1:n_obs), size = 1, replace = TRUE, prob = observation_prob[[a]][s_end, ])
# get the reward
r <- reward[s_start,a]
r_trials[t,h] <- r
r_total <- r_total + discount*r
# update state, belief and discount
s_start <- s_end
b_update <- belief_update(b, model = model_solved$model, act_index = a, obs_index = obs)
b <- b_update
discount <- discount*gamma
}
r_total_trials[t] <- r_total
}
return(list(r_total_trials = r_total_trials,
r_trials = r_trials,
action_trials = action_trials,
state_trials = state_trials,
obs_trials = obs_trials))
}
# Simulate_policy(model_solved = model_POMDPs_solved[[4]],
# policy = "fix",
# n_trial = 1000,
# n_horizon = 60)5.3 Prob_obs_function function
This one is used to calculate \(P(O|belief, a)\).
# compute the P(O|belief, a) function
# Prob_obs_function <- function(b_0, model, act_index){
# n_observations <- length(model$observations) # number of observations
# n_state <- length(model$states) # number of states
#
# transition_prob_act <- model$transition_prob[[act_index]]%>%as.matrix() # transition prob of action
# observation_prob_act <- model$observation_prob[[act_index]]%>%as.matrix() #observation prob of action
# if (observation_prob_act[1,1] == "uniform"){
# observation_prob_act <- matrix(1/n_observations, nrow = n_state, ncol = n_observations)
# }
# # compute the P(O|belief, a) (normalization term)
# prob_obs <- rep(0, n_obs)
# for (o in 1:n_observations){
# sum_1<-0
# sum_2<-0
# for (s_end in 1:n_state){
# for (s_start in 1:n_state){
# x<-b_0[s_start]*transition_prob_act[s_start, s_end]
# sum_1<-sum_1+x
# }
# y<-sum_1*observation_prob_act[s_end, o]
# sum_2<-sum_2+y
# }
# prob_obs[o]<-sum_2
# }
# return(prob_obs)
# }
#
# # Prob_obs_function(b_0 = c(0.25, 0.25, 0.25, 0.25),
# # model= model_POMDPs[[1]],
# # act_index = 4)5.4 Simulate_policy2 function
The main idea in this Simulate_policy2 function is to randomly generate \(state, observation\) at each time step, then the reward of each time step is \(\sum_s R(s,a)*belief(s)\).
# Simulate_policy2 <- function(model_solved, n_trial = 1000, n_horizon = 60,
# policy = "fix", b_0 = "uniform"){
#
# states <- model_solved$model$states
# actions <- model_solved$model$actions
# observations <- model_solved$model$observations
# gamma <- model_solved$model$discount
#
# n_state <- length(states)
# n_obs <- length(observations)
#
# # policy
# if (policy == "optimal"){
# alpha <- model_solved$solution$alpha
# action_index_alpha <- model_solved$solution$action_index_alpha
# }
#
# # initial belief
# if(b_0 == "uniform"){
# b_0 <- rep(1/n_state, n_state)
# }
#
# # reward matrix
# reward <- model_solved$model$reward
# reward <- reward %>%
# mutate(start.state = fct_relevel(start.state, states))%>%
# spread(key = "action", value = "reward")%>%
# select(actions)
#
# # observations matrix
# observation_prob <- model_solved$model$observation_prob
#
# # transition matrix
# transition_prob <- model_solved$model$transition_prob
#
# # ==========Simulation===============
# # run n_trial to get expected rewards
# r_total_trials <- rep(0, n_trial)
# for (t in 1:n_trial){
# discount <- 1
# r_total <- 0
# b <- b_0
# s_start <- sample(c(1:n_state), size = 1, replace = TRUE, prob = b_0) #random generate state: s_start
# for (h in 1:n_horizon){
# # get action index: a
# a <- Action_map(belief = b, policy = policy, alpha, action_index_alpha)
#
# # random generate next state: s_end
# s_end <- sample(c(1:n_state), size = 1, replace = TRUE, prob = transition_prob[[a]][s_start, ])
#
# # random generate observation: obs
# obs <- sample(c(1:n_obs), size = 1, replace = TRUE, prob = observation_prob[[a]][s_end, ])
#
# # get the reward
# r <- b %*% reward[,a]
# r_total <- r_total + discount*r
#
# # update state, belief and discount
# b_update <- belief_update(b, model = model_solved$model, act_index = a, obs_index = obs)
# b <- b_update
# discount <- discount*gamma
# }
#
# r_total_trials[t] <- r_total
# }
# return(r_total_trials)
# }6 ASMBI: Multi-POMDPs
6.1 Step 1: Define and Run POMDPs for each type of hosts.
The defined POMDPs are stored in model_POMDPs and the solved results are stored in model_POMDPs_solved.
# Define hosts features
nOS<-length(transition_prob) #number of OS types
nCrit<-length(reward_list$Linux) #number of critical server type (TRUE and FALSE)
nAdmin<-length(reward_list$Linux$crit_FALSE) #number of admin mode (TRUE and FALSE)
n_type_Hosts <- nOS*nCrit*nAdmin #number of types of hosts
# Define states, observations, actions, and discount
discount<-0.99
states<-c("low" , "medium" , "high", "critical")
actions<-c("auto-patching" , "research-accept", "research-compensate", "remediation")
observations<-c("none", "low", "medium","high","critical")
# Define scenarios when considering errors in observation matrix
errors_scan <- seq(from = 0, to = 0.25, by = 0.05)
errors_name <- paste("error", errors_scan, sep = "_")
n_error <- length(errors_scan) #number of scenarios of observation matrix
# Define start belief, horizons, method, points type for `solve_POMDP2()` function
start<-"uniform"
horizon <- NULL
method <-"grid" #method could be enum, twopass, linsup, witness, incprune, grid, mcgs
params <-list(fg_points = 10000,
fg_type = "initial", # fg_type for grid method and it could be simplex, pairwise, search, initial
epsilon = 1e-1) #set the stop tollerence Run POMDPs: belief points are using “search” method to generate.
# Construct POMDPs and Solve them
## initilize
model_POMDPs <- rep(list(NULL), n_type_Hosts*n_error)
model_POMDPs_solved <- rep(list(NULL), n_type_Hosts*n_error)
model_index <- 0
## loop
pb = txtProgressBar(min = 1, max = (n_type_Hosts*n_error), style = 3) #process bar
for (err in 1:n_error){
for (o in 1:nOS){
for (c in 1:nCrit){
for (a in 1:nAdmin){
#obtain the scenario of errors
Err_name <- errors_name[err]
#obtain hosts' features
OS <- names(reward_list)[o]
Crit <- names(reward_list[[o]])[c]
Admin <- names(reward_list[[o]][[c]])[a]
#obtain the corresponding transitions, observations and rewards
trans_prob <- transition_prob[[OS]]
obs_prob <- observation_prob[[OS]][[Crit]]
reward <- reward_list[[OS]][[Crit]][[Admin]]
#consider errors in observation matrix
obs_prob <- Consider_Error_Observation(obs_matrix = obs_prob,
error_scan = errors_scan[err])
#construct POMDP
model <- pomdp::POMDP(
name = paste(OS, Crit, Admin, Err_name,sep = "_"),
discount = discount,
states = states,
actions = actions,
observations = observations,
start = start,
transition_prob = trans_prob,
observation_prob =obs_prob,
reward = reward)
model_index <- model_index + 1
model_POMDPs[[model_index]]<-model
#progress bar
Sys.sleep(0.02)
setTxtProgressBar(pb, model_index)
#solve POMDP
model_solved <- solve_POMDP2(model = model,
horizon = horizon,
method = method,
verbose = FALSE,
parameter = params)
model_POMDPs_solved[[model_index]] <- model_solved
}
}
}
}Run POMDPs: incremental pruning method
# # Define start belief, horizons, method, points type
# start<-"uniform"
# horizon <- NULL
# method <-"incprune" #method could be enum, twopass, linsup, witness, incprune, grid, mcgs
# params <-list(epsilon = 1e-1) #set the stop tollerence
#
# # Construct POMDPs and Solve them
# ## initilize
# model_POMDPs_2 <- rep(list(NULL), n_type_Hosts*n_error)
# model_POMDPs_solved_2 <- rep(list(NULL), n_type_Hosts*n_error)
# model_index <- 0
#
# ## loop
# pb = txtProgressBar(min = 1, max = (n_type_Hosts*n_error), style = 3) #process bar
#
# for (err in 1:n_error){
# for (o in 1:nOS){
# for (c in 1:nCrit){
# for (a in 1:nAdmin){
#
# #obtain the scenario of errors
# Err_name <- errors_name[err]
# #obtain hosts' features
# OS <- names(reward_list)[o]
# Crit <- names(reward_list[[o]])[c]
# Admin <- names(reward_list[[o]][[c]])[a]
#
# #obtain the corresponding transitions, observations and rewards
# trans_prob <- transition_prob[[OS]]
# obs_prob <- observation_prob[[OS]][[Crit]]
# reward <- reward_list[[OS]][[Crit]][[Admin]]
#
# #consider errors in observation matrix
# obs_prob <- Consider_Error_Observation(obs_matrix = obs_prob,
# error_scan = errors_scan[err])
#
# #construct POMDP
# model <- pomdp::POMDP(
# name = paste(OS, Crit, Admin, Err_name,sep = "_"),
# discount = discount,
# states = states,
# actions = actions,
# observations = observations,
# start = start,
# transition_prob = trans_prob,
# observation_prob =obs_prob,
# reward = reward)
# model_index <- model_index + 1
# model_POMDPs[[model_index]]<-model
# #progress bar
# Sys.sleep(0.02)
# setTxtProgressBar(pb, model_index)
# #solve POMDP
# model_solved_2 <- solve_POMDP2(model = model,
# horizon = horizon,
# method = method,
# verbose = FALSE,
# parameter = params)
# model_POMDPs_solved_2[[model_index]] <- model_solved
# }
# }
# }
# }6.2 Step 2: Shape the outputs
- Get the PG when considering 5 scenarios:
belief_scenarios <- matrix(c(0.25, 0.25, 0.25, 0.25,
1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1),
nrow = 5, ncol = 4, byrow = TRUE)
model_POMDPs_PGs <- lapply(model_POMDPs_solved, function(model_solved){
model_name <- model_solved$model$name
PG <- create_PG2(model = model_solved$model,
alpha = model_solved$solution$alpha,
action_index_alpha = model_solved$solution$action_index_alpha,
belief_matrix = belief_scenarios,
get_action_name = FALSE)
return(list(name = model_name,
PG = PG))
})## Warnings in Other_crit_FALSE_admin_FALSE_error_0 : The belief update is NAs because it is impossible to observe low observation under action auto-patching when the start belief is ( 0,0,1,0 )Warnings in Other_crit_FALSE_admin_TRUE_error_0 : The belief update is NAs because it is impossible to observe low observation under action auto-patching when the start belief is ( 0,0,1,0 )Warnings in Windows_crit_FALSE_admin_FALSE_error_0 : The belief update is NAs because it is impossible to observe low observation under action auto-patching when the start belief is ( 0,0,1,0 )Warnings in Windows_crit_FALSE_admin_TRUE_error_0 : The belief update is NAs because it is impossible to observe low observation under action auto-patching when the start belief is ( 0,0,1,0 )Warnings in Windows_crit_FALSE_admin_TRUE_error_0 : The belief update is NAs because it is impossible to observe low observation under action auto-patching when the start belief is ( 0,0,0,1 )Warnings in Windows_crit_FALSE_admin_TRUE_error_0 : The belief update is NAs because it is impossible to observe medium observation under action auto-patching when the start belief is ( 0,0,0,1 )Warnings in Windows_crit_TRUE_admin_TRUE_error_0 : The belief update is NAs because it is impossible to observe low observation under action auto-patching when the start belief is ( 0,0,1,0 )# assign name to the above list
model_names <-model_POMDPs_PGs%>%map_chr("name")
names(model_POMDPs_PGs) <- model_names
#convert to dataframe
model_POMDPs_PGs_df <- model_POMDPs_PGs%>%{
tibble(
name = map_chr(., "name"),
PG = map(., "PG")
)
}%>%unnest()
# split the name column
name_split <- str_split(model_POMDPs_PGs_df$name, pattern = "_", simplify = T)
model_POMDPs_PGs_df <- model_POMDPs_PGs_df%>%
dplyr::mutate(OS = name_split[, 1],
Criticality = name_split[, 3],
Admin = name_split[, 5],
Error = name_split[, 7])%>%
dplyr::mutate_at(c("Criticality", "Admin"), as.logical)%>%
dplyr::mutate_at("Error", as.numeric)%>%
dplyr::mutate(OS = fct_relevel(OS, c("Linux", "Windows", "Other")))
model_POMDPs_PGs_df%>%filter(OS == "Windows" & Criticality == TRUE & Admin == TRUE)- Get the alpha vectors
model_POMDPs_alpha <- lapply(model_POMDPs_solved, function(model_solved){
alpha_vectors <- model_solved$solution$alpha%>%
dplyr::mutate(action = model_solved$solution$action_index_alpha)%>%
select(action, everything())
model_name <- model_solved$model$name
return(list(name = model_name,
alpha = alpha_vectors))
})
#convert to dataframe
model_POMDPs_alpha_df <- model_POMDPs_alpha%>%{
tibble(
name = map_chr(., "name"),
alpha = map(., "alpha")
)
}%>%unnest()
# split the name column
name_split2 <- str_split(model_POMDPs_alpha_df$name, pattern = "_", simplify = T)
model_POMDPs_alpha_df <- model_POMDPs_alpha_df%>%
dplyr::mutate(OS = name_split2[, 1],
Criticality = name_split2[, 3],
Admin = name_split2[, 5],
Error = name_split2[, 7])%>%
dplyr::mutate_at(c("Criticality", "Admin"), as.logical)%>%
dplyr::mutate_at("Error", as.numeric)%>%
dplyr::mutate(OS = fct_relevel(OS, c("Linux", "Windows", "Other")))
# write in csv
model_POMDPs_alpha_print <- model_POMDPs_alpha_df%>%
mutate(`coeffecient 1` = `coeffecient 1`*-1,
`coeffecient 2` = `coeffecient 2`*-1,
`coeffecient 3` = `coeffecient 3`*-1,
`coeffecient 4` = `coeffecient 4`*-1)
#write.csv(model_POMDPs_alpha_print, file = "Alpha vectors for all POMDPs.csv", row.names = F)
model_POMDPs_alpha_df%>%
filter(OS == "Other" & Criticality == TRUE & Admin == FALSE)- Get the plot of policy graph
7 ASMBI-Simulation: Fixed Policy vs Optimal Policy
if observe nothing, taking action 5(ignore); if observe low and medium, taking action 1; if observe high and critical, taking action 2.
Here, we perform simulations for optimal, fix (1122), OSOM, random and epsilon greedy policeies. The number of horizons is 60. For each policy, the number of simulations is 1000. (The simulations are running for 60 horizons repeated 1000 times)
7.1 Simulation
# initialization
n_POMDPs <- length(model_POMDPs_solved)
policy_compare <- c("optimal", "fix", "OSOM","random","e_greedy")
n_policy <- length(policy_compare)
model_POMDPs_Simulate <- rep(list(NULL), n_POMDPs*n_policy)
loop_index <- 0
## loop
pb = txtProgressBar(min = 1, max = (n_POMDPs*n_policy), style = 3) #process bar
for (n in 1:n_POMDPs){
for (p in 1:n_policy){
#progress bar
loop_index <- loop_index + 1
Sys.sleep(0.01)
setTxtProgressBar(pb, loop_index)
#simulation
policy <- policy_compare[p]
model_name <- model_POMDPs_solved[[n]]$model$name
reward_simulate <- Simulate_policy(model_solved = model_POMDPs_solved[[n]],
policy = policy,
n_trial = 1000,
n_horizon = 60,
b_0 = "random",
seed = 123)
#save results
model_POMDPs_Simulate[[loop_index]] <- list(name = model_name,
policy = policy,
reward = reward_simulate$r_total_trials,
states = reward_simulate$state_trials,
actions = reward_simulate$action_trials)
}
}##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 6%
|
|==== | 7%
|
|===== | 7%
|
|===== | 8%
|
|====== | 9%
|
|====== | 10%
|
|======= | 10%
|
|======= | 11%
|
|======== | 12%
|
|======== | 13%
|
|========= | 13%
|
|========= | 14%
|
|========== | 15%
|
|========== | 16%
|
|=========== | 16%
|
|=========== | 17%
|
|=========== | 18%
|
|============ | 18%
|
|============ | 19%
|
|============= | 19%
|
|============= | 20%
|
|============= | 21%
|
|============== | 21%
|
|============== | 22%
|
|=============== | 23%
|
|=============== | 24%
|
|================ | 24%
|
|================ | 25%
|
|================= | 26%
|
|================= | 27%
|
|================== | 27%
|
|================== | 28%
|
|=================== | 29%
|
|=================== | 30%
|
|==================== | 30%
|
|==================== | 31%
|
|===================== | 32%
|
|===================== | 33%
|
|====================== | 33%
|
|====================== | 34%
|
|====================== | 35%
|
|======================= | 35%
|
|======================= | 36%
|
|======================== | 36%
|
|======================== | 37%
|
|======================== | 38%
|
|========================= | 38%
|
|========================= | 39%
|
|========================== | 39%
|
|========================== | 40%
|
|========================== | 41%
|
|=========================== | 41%
|
|=========================== | 42%
|
|============================ | 42%
|
|============================ | 43%
|
|============================ | 44%
|
|============================= | 44%
|
|============================= | 45%
|
|============================== | 45%
|
|============================== | 46%
|
|============================== | 47%
|
|=============================== | 47%
|
|=============================== | 48%
|
|================================ | 48%
|
|================================ | 49%
|
|================================ | 50%
|
|================================= | 50%
|
|================================= | 51%
|
|================================= | 52%
|
|================================== | 52%
|
|================================== | 53%
|
|=================================== | 53%
|
|=================================== | 54%
|
|=================================== | 55%
|
|==================================== | 55%
|
|==================================== | 56%
|
|===================================== | 56%
|
|===================================== | 57%
|
|===================================== | 58%
|
|====================================== | 58%
|
|====================================== | 59%
|
|======================================= | 59%
|
|======================================= | 60%
|
|======================================= | 61%
|
|======================================== | 61%
|
|======================================== | 62%
|
|========================================= | 62%
|
|========================================= | 63%
|
|========================================= | 64%
|
|========================================== | 64%
|
|========================================== | 65%
|
|=========================================== | 65%
|
|=========================================== | 66%
|
|=========================================== | 67%
|
|============================================ | 67%
|
|============================================ | 68%
|
|============================================= | 69%
|
|============================================= | 70%
|
|============================================== | 70%
|
|============================================== | 71%
|
|=============================================== | 72%
|
|=============================================== | 73%
|
|================================================ | 73%
|
|================================================ | 74%
|
|================================================= | 75%
|
|================================================= | 76%
|
|================================================== | 76%
|
|================================================== | 77%
|
|=================================================== | 78%
|
|=================================================== | 79%
|
|==================================================== | 79%
|
|==================================================== | 80%
|
|==================================================== | 81%
|
|===================================================== | 81%
|
|===================================================== | 82%
|
|====================================================== | 82%
|
|====================================================== | 83%
|
|====================================================== | 84%
|
|======================================================= | 84%
|
|======================================================= | 85%
|
|======================================================== | 86%
|
|======================================================== | 87%
|
|========================================================= | 87%
|
|========================================================= | 88%
|
|========================================================== | 89%
|
|========================================================== | 90%
|
|=========================================================== | 90%
|
|=========================================================== | 91%
|
|============================================================ | 92%
|
|============================================================ | 93%
|
|============================================================= | 93%
|
|============================================================= | 94%
|
|============================================================== | 95%
|
|============================================================== | 96%
|
|=============================================================== | 96%
|
|=============================================================== | 97%
|
|================================================================ | 98%
|
|================================================================ | 99%
|
|=================================================================| 99%
|
|=================================================================| 100%7.2 Shape the results
1. dataframe for rewards
#convert to dataframe
model_POMDPs_Simulate_df <- model_POMDPs_Simulate%>%{
tibble(
name = map_chr(., "name"),
policy = map_chr(., "policy"),
reward = map(., "reward")
)
}%>%unnest()
# split the name column
name_split3 <- str_split(model_POMDPs_Simulate_df$name, pattern = "_", simplify = T)
# dataframe
model_POMDPs_Simulate_df <- model_POMDPs_Simulate_df%>%
dplyr::mutate(OS = name_split3[, 1],
Criticality = name_split3[, 3],
Admin = name_split3[, 5],
Error = name_split3[, 7])%>%
dplyr::mutate_at(c("Criticality", "Admin"), as.logical)%>%
dplyr::mutate_at("Error", as.numeric)%>%
dplyr::mutate(policy = fct_relevel(policy, c("optimal", "fix","OSOM","e_greedy","random")),
OS = fct_relevel(OS, c("Linux", "Windows", "Other")))
model_POMDPs_Simulate_df2. dataframe for average rewards
# take the average
model_POMDPs_Simulate_avg <- model_POMDPs_Simulate_df%>%
group_by(name,OS, Criticality, Admin, Error, policy)%>%
summarise_at("reward", mean)%>%
ungroup()
model_POMDPs_Simulate_avg%>%
filter(policy == "optimal" | policy == "fix")%>%
spread(key = policy, value = reward)%>%
dplyr::mutate(check = abs(optimal) - abs(fix))%>%
filter(check > 0 & Error == 0)3. dataframe for evolutions of states and actions
model_POMDPs_Sim_Evolutions <- model_POMDPs_Simulate%>%{
tibble(
name = map_chr(., "name"),
policy = map_chr(., "policy"),
states = map(., "states"),
actions = map(., "actions")
)
}
# split the name column
name_split4 <- str_split(model_POMDPs_Sim_Evolutions$name, pattern = "_", simplify = T)
model_POMDPs_Sim_Evolutions <- model_POMDPs_Sim_Evolutions%>%
dplyr::mutate(OS = name_split4[, 1],
Criticality = name_split4[, 3],
Admin = name_split4[, 5],
Error = name_split4[, 7])%>%
dplyr::mutate_at(c("Criticality", "Admin"), as.logical)%>%
dplyr::mutate_at("Error", as.numeric)%>%
dplyr::mutate(policy = fct_relevel(policy, c("optimal", "fix","OSOM","e_greedy","random")),
OS = fct_relevel(OS, c("Linux", "Windows", "Other")))
model_POMDPs_Sim_Evolutions7.3 Plots
7.3.1 Costs
- For specific hosts under optimal policy, plot the costs boxplot across observation errors
Admin_names<-c(`TRUE`="Admin Privilege",
`FALSE`="Management")
Critical_names<-c(`TRUE`="Restricted Data",
`FALSE`="Normal Data")
model_POMDPs_Simulate_df%>%
filter(policy == "optimal")%>%
#filter(OS == "Linux" & Criticality == TRUE & Admin == TRUE & policy == "optimal")%>%
ggplot(aes(as.factor(Error), -reward))+
geom_boxplot(aes(fill = OS),
outlier.alpha = 0.5,
outlier.size = 0.5,
varwidth = TRUE)+
facet_grid(Admin ~ OS + Criticality, scales = "free",
labeller = labeller(Admin = as_labeller(Admin_names),
Criticality = as_labeller(Critical_names)))+
scale_fill_nejm()+
theme_bw()+
theme(axis.text.x = element_text(size = 8),
legend.position = "none")+
labs(x = "Error rate of scanning",
y = "5-year costs per host (1000 trials)")
- When error = 0, plot the costs boxplot over different policy
model_POMDPs_Simulate_df%>%
filter(Error == 0)%>%
ggplot(aes(policy, -reward))+
geom_boxplot(aes(fill = OS),
outlier.alpha = 0.8,
outlier.size = 0.5,
varwidth = TRUE)+
facet_grid(Admin ~ OS + Criticality, scales = "free",
labeller = labeller(Admin = as_labeller(Admin_names),
Criticality = as_labeller(Critical_names)))+
scale_fill_nejm()+
theme_bw()+
theme(axis.text.x = element_text(angle=45, vjust=0.6, size =8),
legend.position = "none")+
labs(y = "5-year costs per host (1000 trials)",
x = "Policy")
Exclue some results from the above plot:
# exclue policy = random
model_POMDPs_Simulate_df%>%
filter(Error == 0 & policy != "random")%>%
ggplot(aes(policy, -reward))+
geom_boxplot(aes(fill = OS),
outlier.alpha = 0.5,
outlier.size = 0.5,
varwidth = TRUE)+
facet_grid(Admin ~ OS + Criticality, scales = "free",
labeller = labeller(Admin = as_labeller(Admin_names),
Criticality = as_labeller(Critical_names)))+
scale_fill_nejm()+
theme_bw()+
theme(axis.text.x = element_text(angle=45, vjust=0.6, size =8),
legend.position = "none")+
labs(y = "5-year costs per host (1000 trials)",
x = "Policy")
#ggsave("Costs_boxplot_diff_policy(no_random).pdf", width = 10, height = 6)
# exclude policy = random and fix
model_POMDPs_Simulate_df%>%
filter(Error == 0 & policy != "random" & policy != "fix")%>%
ggplot(aes(policy, -reward))+
geom_boxplot(aes(fill = OS),
outlier.alpha = 0.5,
outlier.size = 0.5,
varwidth = FALSE)+
facet_grid(Admin ~ OS + Criticality, scales = "free",
labeller = labeller(Admin = as_labeller(Admin_names),
Criticality = as_labeller(Critical_names)))+
scale_fill_nejm()+
theme_bw()+
theme(axis.text.x = element_text(angle=45, vjust=0.6, size =8),
legend.position = "none")+
labs(y = "5-year costs per host (1000 trials)",
x = "Policy")
#ggsave("Costs_boxplot_diff_policy(no_random_and_fix).PNG", width = 10, height = 6)
# exclude policy = random, e_greedy
model_POMDPs_Simulate_df%>%
filter(Error == 0 & policy != "random" & policy != "e_greedy")%>%
ggplot(aes(policy, -reward))+
geom_boxplot(aes(fill = OS),
outlier.alpha = 0.5,
outlier.size = 0.5,
varwidth = FALSE)+
facet_grid(Admin ~ OS + Criticality, scales = "free",
labeller = labeller(Admin = as_labeller(Admin_names),
Criticality = as_labeller(Critical_names)))+
scale_fill_nejm()+
theme_bw()+
theme(axis.text.x = element_text(angle=45, vjust=0.6, size =8),
legend.position = "none")+
labs(y = "5-year costs per host (1000 trials)",
x = "Policy")
- Average costs plot: error = 0% and exclude random policy
model_POMDPs_Simulate_avg%>%
filter(Error == 0 & policy != "random")%>%
ggplot(aes(policy, -reward))+
geom_bar(aes(fill = OS),
stat = "identity",
width = 0.5,
color = "black")+
facet_grid(Admin ~ OS + Criticality, scales = "free",
labeller = labeller(Admin = as_labeller(Admin_names),
Criticality = as_labeller(Critical_names)))+
scale_fill_nejm()+
theme_bw()+
theme(axis.text.x = element_text(angle=45, vjust=0.6, size =8),
legend.position = "none")+
labs(y = "5-year average costs per host from 1000 trials",
x = "Policy")
- Average costs plot for only optimal and OSOM: line plot over errors
model_POMDPs_Simulate_avg%>%
filter(policy == "optimal" | policy == "OSOM")%>%
ggplot(aes(x = Error, y = -reward))+
geom_point(aes(colour = OS), size = 1)+
geom_line(aes(group = interaction(OS, policy),colour = OS,
linetype = policy))+
facet_grid(Admin ~ Criticality, scales = "free",
labeller = labeller(Admin = as_labeller(Admin_names),
Criticality = as_labeller(Critical_names)))+
scale_color_nejm()+
theme_bw()+
theme(legend.position = "bottom")+
labs(x = "Observation errors",
y = "5-year average costs per host from 1000 trials")
- Plot only optimal and OSOM: error = 0%, bar plot
model_POMDPs_Simulate_avg%>%
filter(Error == 0 & policy != "random" & policy != "fix" & policy != "e_greedy")%>%
ggplot(aes(policy, -reward))+
geom_bar(aes(fill = OS),
stat = "identity",
width = 0.5,
color = "black")+
geom_text(aes(label = round(-reward, 1)),
size = 2,
vjust = -0.4)+
facet_grid(Admin ~ OS + Criticality, #scales = "free_y",
labeller = labeller(Admin = as_labeller(Admin_names),
Criticality = as_labeller(Critical_names)))+
scale_fill_nejm()+
theme_bw()+
theme(axis.text.x = element_text(angle=45, vjust=0.6, size =8),
legend.position = "none")+
labs(y = "5-year average costs per host from 1000 trials",
x = "Policy")
- Plot savings proportion: optimal vs OSOM when error = 0%
model_POMDPs_Simulate_avg%>%
spread(key = policy, value = reward)%>%
mutate(savings = abs(OSOM)-abs(optimal))%>%
mutate(savings_prop = savings/abs(OSOM))%>%
mutate(OS_short = case_when(
OS == "Windows" ~ "Win",
OS == "Linux" ~ "Linux",
TRUE ~ "Other"
))%>%
mutate(Crit_short = case_when(
Criticality == TRUE ~ "Restricted",
TRUE ~ "Normal"
))%>%
mutate(Admin_short = case_when(
Admin == TRUE ~ "Admin",
TRUE ~ "Mgt"
))%>%
mutate(Type = paste(OS_short, Crit_short, Admin_short, sep = "_"))%>%
mutate(ordering = -as.numeric(OS) + savings_prop)%>%
filter(Error == 0)%>%
ggplot(aes(x = fct_reorder(Type, ordering, .desc = F),
y = 100*savings_prop,
label = savings_prop))+
geom_bar(aes(fill = OS),
colour = "black",
stat="identity",
width=.5)+
geom_text(aes(label = scales::percent(savings_prop)),
position = position_dodge2(width = 0.8),
angle = 0, hjust = -0.3,
size = 3,
check_overlap = F)+
ylim(0, 100)+
coord_flip()+
scale_fill_nejm()+
theme_bw()+
labs(x = "Host Type",
y = "Proportion of Savings (%)")+
theme(legend.position = "none")
- Plot savings proportion: optimal vs fix
model_POMDPs_Simulate_avg%>%
spread(key = policy, value = reward)%>%
mutate(savings = abs(fix)-abs(optimal))%>%
mutate(savings_prop = savings/abs(fix))%>%
mutate(OS_short = case_when(
OS == "Windows" ~ "Win",
OS == "Linux" ~ "Linux",
TRUE ~ "Other"
))%>%
mutate(Crit_short = case_when(
Criticality == TRUE ~ "Restricted",
TRUE ~ "Normal"
))%>%
mutate(Admin_short = case_when(
Admin == TRUE ~ "Admin",
TRUE ~ "Mgt"
))%>%
mutate(Type = paste(OS_short, Crit_short, Admin_short, sep = "_"))%>%
mutate(ordering = -as.numeric(OS) + savings_prop)%>%
filter(Error == 0)%>%
ggplot(aes(x = fct_reorder(Type, ordering, .desc = F),
y = 100*savings_prop,
label = savings_prop))+
geom_bar(aes(fill = OS),
colour = "black",
stat="identity",
width=.5)+
geom_text(aes(label = scales::percent(savings_prop)),
position = position_dodge2(width = 0.8),
angle = 0, hjust = -0.3,
size = 3,
check_overlap = F)+
ylim(-15, 100)+
coord_flip()+
scale_fill_nejm()+
theme_bw()+
labs(x = "Host Type",
y = "Proportion of Savings (%)")+
theme(legend.position = "none")
7.3.2 Plot the evolution of states
- Here, we create
Plot_evolutionfunction to plot a heat map of evolutions for specific policy and type of hosts.
# plot evolutions for specific policy
Plot_evolution <- function(data_evolution, target = "states", trial_limit = 1000, plot_ratio = 0.1,
OS = "Linux",
Criticality = TRUE,
Admin = TRUE,
Error = 0,
policy = "optimal"){
os <- OS
crit <- Criticality
admin <- Admin
err <- Error
pol <- policy
limit <- trial_limit
# symbolize
target <- as.name(target)
# enquotes
target <- enquo(target)
# string
target_name <- quo_name(target)
data <- data_evolution%>%
filter(OS == os & Criticality == crit & Admin == admin & Error == err & policy == pol)
data_matrix <- data[[target_name]][[1]]
data_df <- setNames(reshape2::melt(data_matrix), c("trial", "horizons", target_name))
plot <- data_df%>%
filter(trial <= limit)%>%
ggplot(aes(x = horizons, y = trial, fill = !!target) )+
geom_tile(colour="white",size=0.01)+
scale_fill_gradient2(low = "white", high = "black")+
scale_y_continuous(expand=c(0,0))+
scale_x_continuous(expand = c(0,0),
breaks = seq(0, 60, by = 6))+
coord_fixed(ratio = plot_ratio)+
theme_bw()
#theme(plot.margin = unit(c(0,0,1,1), "mm"))
return(plot)
}Plot_evolution_all: plot evolutions for type of host
Plot_evolution_all <- function(data_evolution, target = "states", trial_limit = 1000, plot_ratio = 0.1,
OS = "Linux",
Criticality = TRUE,
Admin = TRUE,
Error = 0,
remove_random = TRUE,
remove_fix = TRUE){
os <- OS
crit <- Criticality
admin <- Admin
err <- Error
limit <- trial_limit
# symbolize
target <- as.name(target)
# enquotes
target <- enquo(target)
# string
target_name <- quo_name(target)
data <- data_evolution%>%
filter(OS == os & Criticality == crit & Admin == admin & Error == err)
policy <- data$policy
data_list <- lapply(policy, function(p){
data_p <- data%>%filter(policy == p)
data_matrix <- data_p[[target_name]][[1]]
data_df <- setNames(reshape2::melt(data_matrix), c("trial", "horizons", target_name))
data_df <- data_df%>%mutate(policy = p)
return(data_df)
})
data_all <- data_list%>%map_df(`[`) #convert list to dataframe
if(remove_random & remove_fix){
plot <- data_all%>%
filter(trial <= limit & policy != "random" & policy != "fix")%>%
ggplot(aes(x = horizons, y = trial, fill = !!target) )+
facet_grid(~policy)+
geom_tile(colour="white",size=0.01)+
scale_fill_gradient2(low = "white", high = "black")+
scale_y_continuous(expand=c(0,0))+
scale_x_continuous(expand = c(0,0),
breaks = seq(0, 60, by = 6))+
coord_fixed(ratio = plot_ratio)+
theme_bw()+
theme(axis.text.x = element_text(size = 5),
legend.position = "bottom")
}else if (remove_random & !remove_fix){
plot <- data_all%>%
filter(trial <= limit & policy != "random")%>%
ggplot(aes(x = horizons, y = trial, fill = !!target) )+
facet_grid(~policy)+
geom_tile(colour="white",size=0.01)+
scale_fill_gradient2(low = "white", high = "black")+
scale_y_continuous(expand=c(0,0))+
scale_x_continuous(expand = c(0,0),
breaks = seq(0, 60, by = 6))+
coord_fixed(ratio = plot_ratio)+
theme_bw()+
theme(axis.text.x = element_text(size = 8),
legend.position = "bottom")
}else if(!remove_random & remove_fix){
plot <- data_all%>%
filter(trial <= limit & policy != "fix")%>%
ggplot(aes(x = horizons, y = trial, fill = !!target) )+
facet_grid(~policy)+
geom_tile(colour="white",size=0.01)+
scale_fill_gradient2(low = "white", high = "black")+
scale_y_continuous(expand=c(0,0))+
scale_x_continuous(expand = c(0,0),
breaks = seq(0, 60, by = 6))+
coord_fixed(ratio = plot_ratio)+
theme_bw()+
theme(axis.text.x = element_text(size = 8),
legend.position = "bottom")
}else{
plot <- data_all%>%
filter(trial <= limit)%>%
ggplot(aes(x = horizons, y = trial, fill = !!target) )+
facet_grid(~policy)+
geom_tile(colour="white",size=0.01)+
scale_fill_gradient2(low = "white", high = "black")+
scale_y_continuous(expand=c(0,0))+
scale_x_continuous(expand = c(0,0),
breaks = seq(0, 60, by = 6))+
coord_fixed(ratio = plot_ratio)+
theme_bw()+
theme(axis.text.x = element_text(size = 8),
legend.position = "bottom")
}
return(plot)
}Example:
# Plot_evolution(model_POMDPs_Sim_Evolutions,
# target = "states", trial_limit = 1000, plot_ratio = 0.1,
# OS = "Linux",
# Criticality = TRUE,
# Admin = FALSE,
# Error = 0,
# policy = "optimal")
#
# Plot_evolution(model_POMDPs_Sim_Evolutions,
# target = "states", trial_limit = 1000,plot_ratio = 0.1,
# OS = "Linux",
# Criticality = TRUE,
# Admin = FALSE,
# Error = 0,
# policy = "fix")
#ggsave("test1.PNG", width = 4, height = 5)- Plot evolution of states for all policies: Linux, critical=true, admin=true
Plot_evolution_all(model_POMDPs_Sim_Evolutions,
target = "states", #could be states or actions
trial_limit = 200,
plot_ratio = 0.5,
OS = "Linux", #could be Linux, Windows, Other
Criticality = TRUE,
Admin = TRUE,
Error = 0,
remove_random = T,
remove_fix = F)
#ggsave("Evolution_states_policies(Linux_crit_True_admin_True).PNG", width = 9, height = 5)
#ggsave("Evolution_states_policies(Linux_crit_True_admin_True)_no_random.PNG", width = 9, height = 5)
#ggsave("Evolution_states_policies(Linux_crit_True_admin_True)_fist100_no_random.PNG", width = 9, height = 5)
#ggsave("Evolution_states_policies(Linux_crit_True_admin_True)_fist200_no_random.pdf", width = 9, height = 5)- Plot evolution of actions for all policies: Linux, critical=true, admin=true
Plot_evolution_all(model_POMDPs_Sim_Evolutions,
target = "actions", #could be states or actions
trial_limit = 200,
plot_ratio = 0.5,
OS = "Linux", #could be Linux, Windows, Other
Criticality = TRUE,
Admin = TRUE,
Error = 0,
remove_random = T,
remove_fix = F)
#ggsave("Evolution_states_policies(Linux_crit_True_admin_True).PNG", width = 9, height = 5)
#ggsave("Evolution_states_policies(Linux_crit_True_admin_True)_no_random.PNG", width = 9, height = 5)
#ggsave("Evolution_states_policies(Linux_crit_True_admin_True)_fist100_no_random.PNG", width = 9, height = 5)
#ggsave("Evolution_states_policies(Linux_crit_True_admin_True)_fist200_no_random.PNG", width = 9, height = 5)- Plot evolution of states for all policies: Linux, critical=true, admin=False
Plot_evolution_all(model_POMDPs_Sim_Evolutions,
target = "states", #could be states or actions
trial_limit = 200,
plot_ratio = 0.5,
OS = "Linux", #could be Linux, Windows, Other
Criticality = TRUE,
Admin = FALSE,
Error = 0,
remove_random = T,
remove_fix = F)
#ggsave("Evolution_states_policies(Linux_crit_True_admin_False).PNG", width = 9, height = 5)
#ggsave("Evolution_states_policies(Linux_crit_True_admin_False)_fist200_no_random.eps", width = 9, height = 5)
#ggsave("Evolution_states_policies(Linux_crit_True_admin_False)_no_random.PNG", width = 9, height = 5)- Plot evolution of states for all policies: Windows, critical=true, admin=true
Plot_evolution_all(model_POMDPs_Sim_Evolutions,
target = "states", #could be states or actions
trial_limit = 200,
plot_ratio = 0.5,
OS = "Windows", #could be Linux, Windows, Other
Criticality = TRUE,
Admin = TRUE,
Error = 0,
remove_random = T,
remove_fix = F)
- Plot evolution of actions for all policies: Windows, critical=true, admin=true
Plot_evolution_all(model_POMDPs_Sim_Evolutions,
target = "actions", #could be states or actions
trial_limit = 200,
plot_ratio = 0.5,
OS = "Windows", #could be Linux, Windows, Other
Criticality = TRUE,
Admin = TRUE,
Error = 0,
remove_random = T,
remove_fix = F)
- Compare Windows and Linux when policy = optimal
Plot_evolution(model_POMDPs_Sim_Evolutions,
target = "states",
trial_limit = 200,
plot_ratio = 0.5,
OS = "Linux",
Criticality = TRUE,
Admin = TRUE,
Error = 0,
policy = "optimal")
Plot_evolution(model_POMDPs_Sim_Evolutions,
target = "states",
trial_limit = 200,
plot_ratio = 0.5,
OS = "Windows",
Criticality = TRUE,
Admin = TRUE,
Error = 0,
policy = "optimal")
7.4 Observations (example)
data_imp<-read.csv("Data for POMDP/Alldata for R.csv",header=TRUE,na.strings=c(""))
# Add IP as new column in to the host-based data
IP_list<-read.csv("Data for POMDP/IPaddressList.csv",header=TRUE,na.strings=c(""))
data_imp$IP<-rep(IP_list$IP,21)Look at the state of each host over 21 month
Host_occurance<-data_imp%>%select(c("IP","HostState","ScannedPeriod"))%>%
arrange(IP, ScannedPeriod)%>%
spread(ScannedPeriod, HostState)
#print(Host_occurance, nrow = 15, ncol = 22)
#write.csv(Host_occurance, file = "Host_occurance_over21month.csv")Basically, we’re looking for the proportion of disappeared time in the 21 scanned period for each unique host. Then, we intend to obtain the mean of the proportion of disappeared time for each type of hosts.
For each unique host or IP, generate the vector of occurrence over months. Then, the proportion of missing is the number of 0’s/the totoal number of months. After that, we could get a set of proportions of missing for each OS type of hosts.
Prop_table<-data_imp%>%
dplyr::group_by(IP, HostState)%>%
dplyr::summarise(N = n())%>%
dplyr::mutate(Proportion = N/sum(N))%>%
dplyr::filter(HostState>0)%>%
dplyr::mutate(Occurrence = sum(Proportion), Disappearance = 1-Occurrence)
# generate the hosts IPs that alway exist over 21 months
# Note that these IPs have done the imputation procedure, which means the number of hosts always existing is larger than the original one (8582 vs. 2288)
IP_alwaysExisting<-subset(Prop_table, Occurrence == 1)
IP_alwaysExisting<-IP_alwaysExisting[!duplicated(IP_alwaysExisting$IP), ]
# Go back to the proportions of missing
# merge the "OS.Type", "Critical","GeneralPurpose", "Host.Type", "AdminPrivilege"
data_imp_sub<-subset(data_imp, ScannedPeriod ==1,
select = c("IP", "OS.Type", "Critical","GeneralPurpose", "Host.Type", "AdminPrivilege"))
Prop_table<-merge(Prop_table, data_imp_sub,by = c("IP"))
Prop_missing<-subset(Prop_table,
select = c("IP", "Occurrence", "Disappearance", "OS.Type", "Critical","GeneralPurpose",
"Host.Type", "AdminPrivilege"))
Prop_missing<-Prop_missing[!duplicated(Prop_missing$IP), ]
as.tibble(Prop_missing)Then, we could draw the box plot for each OS type of hosts
# check the distribution of the disappeared percent for different OS and different data type
Prop_missing_means <- aggregate(Disappearance ~ OS.Type + Critical, Prop_missing, mean)
g<-ggplot(Prop_missing, aes(OS.Type, Disappearance))+
geom_boxplot(aes(fill = factor(Critical)), colour = "black",
position=position_dodge(.8))+
stat_summary(fun.y = mean,
colour = "black",
geom = "point",
aes(group = factor(Critical)),
shape = 18,
size = 3,
position=position_dodge(.8)) +
# stat_summary(aes(label = ..y.., group = factor(Critical)),
# fun.y = mean,
# geom = "text")+
geom_text(data = Prop_missing_means,
aes(label = round(Disappearance,3),
group = factor(Critical)),
position = position_dodge(.8),
vjust = -0.9,
size = 3)+
theme_bw()+
#scale_fill_brewer(palette = "Reds")+
labs(x="Operating Systems",
y="Disappearance proportion",
fill = "Restricted")+
theme(legend.position = "bottom")
print(g)
Then, we can calculate the average, median of disappeared proportion and corresponding standard deviation for each group shown above.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.09524 0.42857 0.45308 0.85714 0.95238# For each group :
df<-Prop_missing%>%
dplyr::group_by(OS.Type, Critical)%>%
dplyr::summarise(Average = mean(Disappearance), Median = median(Disappearance),
SD = sd(Disappearance) )
arrange(df, Critical)